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ABSTRACT 

We propose a two-stage algorithm for the super-resolution of 
MR images from their low-frequency k-space samples. In the 
first stage we estimate a resolution-independent mask whose 
zeros represent the edges of the image. This builds off recent 
work extending the theory of sampling signals of finite rate 
of innovation (FRI) to two-dimensional curves. We enable 
its application to MRI by proposing extensions of the signal 
models allowed by FRI theory, and by developing a more ro¬ 
bust and efficient means to determine the edge mask. In the 
second stage of the scheme, we recover the super-resolved 
MR image using the discretized edge mask as an image prior. 
We evaluate our scheme on simulated single-coil MR data ob¬ 
tained from analytical phantoms, and compare against total 
variation reconstructions. Our experiments show improved 
performance in both noiseless and noisy settings. 

Index Terms — Super-resolution, MRI, Finite Rate of In¬ 
novation Curves 

1. INTRODUCTION 

The availability of high-resolution MRI data can greatly facil¬ 
itate early diagnosis by enabling the detection and character¬ 
ization of subtle and clinically significant lesions Q. How¬ 
ever, the recovery of very high-resolution MRI data is often 
challenging, mainly due to the slow nature of MRI acquisi¬ 
tion, subject motion, and the rapid decrease in SNR with res¬ 
olution. For example, it is common practice to acquire low- 
resolution data in MR spectroscopic imaging since the acqui¬ 
sition of higher k-space samples comes with a heavy SNR 
penalty, which is undesirable for imaging metabolites at very 
low concentrations. At the same time, the use of low spa¬ 
tial resolution results in the leakage of strong signals such 
as water and fat to other spatial regions, thus distorting the 
metabolite signals. 

Several super-resolution and off-the-grid methods were 
introduced recently to estimate parametric signals with finite 
number of unknowns from their low-frequency samples. For 
example, the locations and amplitudes of a finite number of 
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Dirac delta functions can be reliably estimated from its low- 
frequency Fourier samples using continuous total variation 
minimization j^, or by estimating an annihilation filter and 
then determining its roots 0- Unfortunately, the represen¬ 
tation using finite number of basis functions is not efficient 
for the representation of piecewise polynomial images. For 
example, the gradient of a piecewise constant image is non¬ 
zero on a curve, which cannot be represented as a finite linear 
combination of Dirac delta functions. 

Recently, Pan et al. introduced an complex analytic sig¬ 
nal model for continuous domain two-dimensional signals, 
whose derivatives are supported on a curve 0. Specifically, 
they assume the curve is the zero level-set of a function band- 
limited to a rectangular region in the Fourier domain. The 
authors showed that the Fourier transform of the curve model 
will annihilate the Fourier coefficients of the signal deriva¬ 
tives. This property enables them to extend the FRI model 0 
to multidimensional signals. 

The main focus of this paper is to extend the multidimen¬ 
sional FRI model 0 to enable super-resolution MRI. While 
the multidimensional FRI model 0 is very powerful, it has 
some limitations that restricts its direct applicability to MRI. 
First of all, the complex analytic signal model introduced in 
0 is too restrictive and is not applicable to most MR im¬ 
ages. We generalize the signal model to piecewise polyno¬ 
mial and harmonic functions, which are better suited to repre¬ 
sent practical signals. The annihilation conditions central to 
the scheme are not exactly satisfied in the presence of model 
mismatch and noise. A Cadzow iterative procedure was used 
in 0 to denoise the data before estimating the curve coef¬ 
ficients. This approach is computationally demanding, and 
requires some knowledge of the underlying model order. We 
introduce a novel algorithm based on averaging vectors in the 
null space of the equations that is robust to noise and model- 
mismatch, and is computationally efficient. 

Our reconstruction algorithm proceeds in two steps: First, 
we estimate a super-resolved spatial mask whose zeros corre¬ 
spond to the edges in the image. Second, we discretize the 
mask at the desired resolution and use the discrete spatial 
weights in a weighted total variation. We compare the effi¬ 
cacy of this two step strategy against classical discrete total 
variation regularized recovery of numerical phantoms from 
their exact Fourier samples. 



2. SUPER-RESOLUTION OF EDGE IMAGES 

We will show under certain conditions it is possible to super¬ 
resolve the edges of an image from its low-frequency Fourier 
samples. In general we cannot hope to recover the edge set 
unless we assume certain constraints on its geometry. As in 
0 we assume the edge sets to be the zero set of a bandlimited 
periodic trigonometric polynomial /i(r) on [0,1]^, 

C : = 0 (1) 
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M(r) 


The above proposition shows that in principle it is possi¬ 
ble to recover the edge set C : {/i = 0} of by solving the 
linear system of equations 0 for /i. This is provided we have 
a sufficient number of low-pass samples to make the system 
0 determined. We investigate methods for solving for the 
filter coefficients yu in the following section. 

We now consider the case where the image is assumed 
to be piecewise linear, meaning it can be written as a linear 
combination of functions of the form 


5(r) 
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where A is any finite index set, and yu[k] are any complex 
coefficients. We call such sets C trigonometric curves. As 
noted in the set of trigonometric curves have a rich topol¬ 
ogy, and for large enough bandwidth, they can approximate 
arbitrary curves to any desired accuracy. 

As a first approximation, we consider images that are 
piecewise constant, meaning the image can be expressed as 
finite linear combinations of functions of the form 


. /1 if r = (a;, y) 

else. 


where is a simple region in [0,1]^ with piecewise smooth 
boundary dVt given as {/i = 0} for some /i as in ([T]). The 
Fourier transform of is given as 
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where the last two equations follow by Green’s theorem. Un¬ 
der the Fourier domain relations 


where L is any affine function L{r) = a^r + b, for a G 
C^, b e C. Intuitively, any second derivative of g should 
vanish except on dfl : {/i(r) = 0}, where it will act like the 
derivative of a Dirac. Accordingly, we can show z/ = /i^ is an 
annihilating mask for any second derivative of g, since both 
a and Vz/ = 2/iV/i vanish on dfl. The following proposition 
expresses this fact in the Fourier domain: 

Proposition 2.2. Let g = L • where d‘^L = 0, for some 
second order differential operator d^. Then the Fourier trans¬ 
form ofh = d^g is annihilated by convolution with /i^[k] on 
its support A', that is 

/^^[k] h{uj — 27rk) = 0, for all u; G 

keA' 

Note that a similar claim holds with the Laplacian A in 
place of d‘^ and when L is any harmonic function, AL = 0. 
Therefore, images that are piecewise harmonic will satisfy the 
same annihilation property in Prop. |2.2| 

3. ALGORITHMS 

3.1. Computation of an annihilating filter 


5a;ln(r)o - jwa;ln(w), dylii{r) o -jwyln(w) 

the partial derivatives of can be interpreted (in a distribu¬ 
tional sense) as a continuous stream of Diracs in the spatial 
domain supported on dQ. Thus, formally, we should have 
F ' = F ' dy^n = 0^ so we say g acts as an an¬ 

nihilating mask for the partial derivatives of 1 q. This can be 
established rigorously in the Fourier domain to give the fol¬ 
lowing: 

Proposition 2.1. Let be as above with boundary dVt given 
as a trigonometric curve C : {/i(r) = 0}. Then the Fourier 
transforms of the partial derivatives ofl^t are annihilated by 
convolution with yu[k], that is 

/i[k] f{uj — 27rk) = 0, for all u; G (4) 

kGA 

where f{u;) = orf{uj) := 


Suppose we wish to find filter coefficients c G supported 
in A that annihilate the low-frequency Fourier coefficients of 
the image derivatives dj G C^, j = 1,..., n, that is 

{dj * c)[k] = 0 for all kGA, j = 1, n. 

This is equivalent to finding a non-zero vector c in the 
nullspace of the each block Toeplitz matrix G 
corresponding to 2-D convolution with d^. Setting T = 
[Tf,..., T^]^, one approach would be to solve the least 
squares problem 

min||Tc|| 2 , subject to c[0] = 1. (6) 

c 

While computationally efficient, this method is highly sensi¬ 
tive to noise. Moreover, when the model order M is larger 
than strictly necessary, this can lead to spurious zeros in the 
reconstructed mask g = [c]. 






In the case of noisy measurements, the authors of 0 sug¬ 
gest first performing a Cadzow denoising Q of the Toeplitz 
matrix T, then solving However, the computational de¬ 
mands of Cadzow’s algorithm are significant in this case, 
since each iteration of the algorithm requires finding the S VD 
of a large convolution matrix. Moreover, we found Cadzow’s 
algorithm was unable to correct for spurious zeros unless the 
model order reduction was chosen correctly, which is difficult 
to know in advance. 

To reduce the impact of noise and avoid spurious zeros 
in our annihilating mask, we propose a procedure that aver¬ 
ages a collection of elements in the (approximate) nullspace 
of T. First, consider the noiseless overdetermined case where 
T is low-rank. To obtain a basis for the nullspace of T we 
compute its SVD, T = U5]V^, and collect the columns 
of V corresponding to zero singular values, labeling them 
From each we can construct an annihilating mask 
= T~^[wi]. Individually each may contain extra ze¬ 
ros outside C, but the intersection of their zero sets is likely 
only to contain C, due to orthogonality constraints. Hence 
we may compute an improved annihilating mask p by form¬ 
ing the sum-of-squares average 
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In the presence of noise T will only be approximately low- 
rank, and in general we may not know the minimal model 
order representing the edge image, meaning we cannot pre¬ 
cisely identify the nullspace vectors {v^}. Therefore we pro¬ 
pose collecting the {v^} corresponding to singular values 
less than some threshold r = 5 ’di, where 6 G (0,1] is a con¬ 
stant depending on the noise level. In this case, forming the 
averaged mask in 0 has the additional benefit of smoothing 
out variations in each mask due to noise. 

In Fig.we compare the proposed method for computing 
an edge mask versus the least-squares and Cadzow denoising 
methods proposed in 0 on the Shepp-Logan phantom in both 
low noise and high noise settings. Here we corrupt the k- 
space samples with mean-zero additive white Gaussian noise 
to have the indicated signal-to-noise ratio (SNR). We observe 
that the proposed method gives significantly improved masks 
in both cases. 


3.2. Reconstruction algorithm 

Following 0, we propose a two-stage algorithm for recov¬ 
ering MR images. The first stage is to estimate an annihilat¬ 
ing mask fi using the procedure in the previous subsection. 
Note that this first step can be carried out in a resolution- 
independent manner—the mask fi records the edge locations 
with infinite resolution. In the second stage, we discretize the 
mask /i to the desired resolution as a set of spatial weights 



(a) Least-Squares (b) Cadzow (c) Proposed 



(d) Least-Squares (e) Cadzow (f) Proposed 


Fig. 1: Edge annihilation masks for Shepp-Logan phantom 
using 61^ k-space samples. Top row shows a low-noise set¬ 
ting (SNR=25 dB), bottom row shows a high-noise setting 
(SNR=20 dB). The proposed method yields a significantly 
improved mask in both settings. 


and perform a weighted total variation (TV) recovery: 

min||Ax-b||2 + A||W^-|Vx|||i (8) 

X 

where x the discrete image to be recovered, A is the Fourier 
undersampling operator, b is the vector of noisy low-pass 
Fourier measurements, A is a regularization parameter, and 
I Vx| denotes the magnitude of the discrete gradient of x. In 
the case where a higher order signal model is assumed, we 
could also replace the weighted TV penalty with a weighted 
higher degree TV (HDTV) penalty 0. 

4. RESULTS 

In this section we demonstrate the utility of our proposed 
recovery scheme for single-coil MR image recovery from 
low-pass k-space samples. For simplicity, we focus on the 
first-order case, using piecewise constant phantoms. In order 
to obtain resolution-independent simulations, we calculate 
k-space data directly using the analytical MRI phantoms 
derived in 0 - We simulate a noisy acquisition by adding 
complex white Gaussian noise with standard deviation a to 
the computed k-space samples. 

We collect Cartesian center k-space samples located 
within a rectangular window [—2A, 2A] x [—21/, 21/], and 
solve for an edge mask using the procedure outlined in Sec. 
US We then recover the image using the weighted TV 
scheme 0, and compare against standard TV, i.e. (0 with 
= 1. In both cases we tune the regularization parame¬ 
ter A to optimize the signal-to-noise ratio (SNR), defined as 









(a) Fully sampled, zoom 



(b) TV, SNR=16.6dB (c) Proposed, SNR=21.3dB 



(d) TV, SNR=14.6dB (e) Proposed, SNR=20.1dB 



(f) Fully sampled, zoom 
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(g) TV, SNR=19.1dB (h) Proposed, SNR=19.0dB 



(i) TV, SNR=16.9dB (j) Proposed, SNR=17.1dB 


Fig. 2 : (Top row) Shepp-Logan phantom reconstructed onto a 256 x 256 spatial grid from 65 x 49 center k-space samples 20-fold 
undersampling) with no added noise in (b)&(c) and 25dB added noise (d)&(e). (Bottom row) Brain phantom reconstructed onto 256x256 
spatial grid reconstructed from 97 x 97 center k-space samples 7-fold undersampling) with no added noise in (g)&(h) and 30dB added 
noise (i)&(j). Note the proposed scheme gives sharper reconstructions with fewer aliasing artifacts over standard TV. 


SNR = 201ogio(||xo||2/||x — X 0 II 2 ) where x is the recon¬ 
structed image, and the ground truth xq is a TV reconstruction 
of the image from fully-sampled k-space. 

Figure shows the reconstructions obtained at a resolu¬ 
tion of 256 X 256 from 65 x 49 = 3185 k-space samples in the 
case of the Shepp-Logan phantom, and from 97 x 97 = 9409 
noisy k-space samples for the brain phantom, in both noise¬ 
less and noisy setting. We observe that in the standard TV 
reconstructions, using a low regularization weight recovers 
much of the fine detail but results in significant aliasing arti¬ 
facts. However, increasing the regularization weight enough 
to eliminate the aliasing results in oversmoothing. In contrast, 
performing a weighted TV with an edge mask allows us to si¬ 
multaneously recover the fine details while also suppressing 
aliasing artifacts. 

5. CONCLUSIONS 

We extend the theory in Q to higher-order signal models, 
and propose a more robust and efficient means of computing 
super-resolved edge images. Our experiments on analytical 
phantoms demonstrate the potential of the proposed scheme 
for super-resolution MRI. 
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